/* --------------------------------------
SUMMARY: 	- Estimates Discrete Choice Models
			- Exports Results to Table
AUTHOR: Thimo De Schouwer (KU LEUVEN)
STATA VERSION: SE 18.5
--------------------------------------- */

set more off, perm

clear all
drop _all

** Recover directory
cd "${dir}\DataCode\Data\Constructed\LISS"
import delimited Clean_LISS.txt, clear

*=========================*
* 1. Reshaping data 	  *
*=========================*
* Job A characterized by: whA flA tcA mwA mwageA
* Job B characterization differs by experiment
* 2.1 First reshape: one observation per experiment

** JOB A (based on actual job)
forvalues i = 1/10{
	gen tca`i' = telecommutea
	gen wha`i' = workhoursa
	gen fla`i' = flexibilitya
	gen mwa`i' = meaninga
	gen hwagea`i' = hrlywagea
	gen mwagea`i' = monthlywagea
}

** JOB B (differs by experiment)
* 1. Workhours
foreach i in 2 3 4 8 9 10{
	gen whb`i' = workhoursa
}

* 2. Flexibility
foreach i in 1 3 4 5 7 10{
	gen flb`i' = flexibilitya
}

* 3. Meaning
foreach i in 4 5 8 10{
	gen mwb`i' = meaningb
}

foreach i in 1 2 3 6 7 9{
	gen mwb`i' = meaninga
}

* 4. Telecommute
foreach i in 3 7 9 10{
	gen tcb`i' = telecommuteb
}

foreach i in 1 2 4 5 6 8 {
	gen tcb`i' = telecommutea
}

* 5. Monthly Wages & Hourly Wages
forvalues i = 1/10{
	gen mwageb`i' = wage_b`i'
}

forvalues i = 1/10{
	gen hwageb`i' = mwageb`i' / (whb`i'*(52/12))
}


* 6. Choices
forvalues i = 1/10{
	gen choicea`i' = 1 if hce`i' == 1
		replace choicea`i' =  0 if hce`i' == 2
	gen choiceb`i' = 1 if hce`i' == 2
		replace choiceb`i' =  0 if hce`i' == 1
}


reshape long wha whb mwagea mwageb hwagea hwageb mwa mwb tca tcb fla flb choicea choiceb, i(nomem_encr) j(experiment)

* 1.2 Second reshape: one observation per choice
reshape long wh mw tc fl choice hwage mwage, i(nomem_encr experiment) j(choicealt a b) string
sort nomem_encr experiment

* Generating log wages and earnings
gen ln_mwage = ln(mwage)
gen ln_hwage = ln(hwage)

* Binary flexibility and hour categories
recode wh (38 = 0) (32 = 1) (20 = 2)
recode fl (0 = 0) (1 2 3 = 1)


*=======================================*
* 2. 	Construct Interaction Terms		*
*=======================================*
* Interaction terms: Job Attribute x Young Child
rename children child 
gen child_wh1 = (child * wh) if wh == 1
	replace child_wh1 = 0 if child_wh1 == .
	
gen child_wh2 = (child * wh) if wh == 2
	replace child_wh2 = 0 if child_wh2 == .
	
gen child_fl = fl * child
gen child_tc = tc * child
gen child_mw = mw * child


label var child_wh1 "Child $\times$ Long PT"
label var child_wh2 "Child $\times$ Short PT"
label var child_tc "Child $\times$ Telecommute"
label var child_mw "Child $\times$ Meaning"
label var child_fl "Child $\times$ Schedule Freedom"


global interactions child_wh1 child_wh2 child_tc child_mw child_fl

*===========================================================*
* 3. 	Estimating Discrete Choice Models 					*
*		Preferences Vary over Wage Distribution			  	*
*===========================================================*
drop if choice == .		// drop empty experiments that weren't conducted

* Declare choice method data
cmset nomem_encr experiment choicealt

* Mean wages by gender
forvalues g = 0/1{
	mean mwage if woman == `g'
	gl meanwage_`g' = r(table)[1,1]
}

eststo coefs: cmclogit choice ib0.fl ib0.tc ib0.mw ln_hwage i.wh ib0.mw#ib0.tc ib0.mw#ib0.fl $interactions if woman == 0
eststo coefs: cmclogit choice ib0.fl ib0.tc ib0.mw ln_hwage i.wh ib0.mw#ib0.tc ib0.mw#ib0.fl $interactions if woman == 1


* Estimating the model (NB: equivalent to not doing second recast and estimating logit choiceA on difference in amenities between jobs)
forvalues g = 0/1{
	eststo coefs_`g': cmclogit choice fl tc mw ln_hwage i.wh $interactions if woman == `g', vce(cluster nomem_encr)
	
		* WTP in Euro: no children
		eststo EUR_WTP_SE_nc_`g': nlcom (EUR_WTP_SE_nc_fl_`g': (1 - exp(- (_b[fl]) / _b[ln_hwage])) * ${meanwage_`g'}) ///
				(EUR_WTP_SE_nc_tc_`g': (1 - exp(- (_b[tc] / _b[ln_hwage]))) * ${meanwage_`g'}) ///
				(EUR_WTP_SE_nc_mw_`g': (1 - exp(- (_b[mw] / _b[ln_hwage]))) * ${meanwage_`g'}) ///
				(EUR_WTP_SE_nc_wh1_`g': (1 - exp(- (_b[1.wh] / _b[ln_hwage]))) * ${meanwage_`g'}) ///
				(EUR_WTP_SE_nc_wh2_`g': (1 - exp(- (_b[2.wh] / _b[ln_hwage]))) * ${meanwage_`g'}), post
		
		est restore coefs_`g'
		* WTP in Euro: children
		eststo EUR_WTP_SE_c_`g': nlcom (EUR_WTP_SE_c_fl_`g': (1 - exp(- (_b[fl] + _b[child_fl]) / _b[ln_hwage])) * ${meanwage_`g'}) ///
				(EUR_WTP_SE_c_tc_`g': (1 - exp(- ((_b[tc] + _b[child_tc]) / _b[ln_hwage]) )) * ${meanwage_`g'}) ///
				(EUR_WTP_SE_c_mw_`g': (1 - exp(- ((_b[mw] + _b[child_mw]) / _b[ln_hwage]) )) * ${meanwage_`g'}) ///
				(EUR_WTP_SE_c_wh1_`g': (1 - exp(- ((_b[1.wh] + _b[child_wh1]) / _b[ln_hwage]) )) * ${meanwage_`g'}) ///
				(EUR_WTP_SE_c_wh2_`g': (1 - exp(- ((_b[2.wh] + _b[child_wh2]) / _b[ln_hwage]) )) * ${meanwage_`g'}), post
		
		est restore coefs_`g'		
		* WTP in Euro: no children
		eststo WTP_SE_nc_`g': nlcom (EUR_WTP_SE_nc_fl_`g': (1 - exp(- (_b[fl]) / _b[ln_hwage]))) ///
				(WTP_SE_nc_tc_`g': (1 - exp(- (_b[tc] / _b[ln_hwage])))) ///
				(WTP_SE_nc_mw_`g': (1 - exp(- (_b[mw] / _b[ln_hwage])))) ///
				(WTP_SE_nc_wh1_`g': (1 - exp(- (_b[1.wh] / _b[ln_hwage])))) ///
				(WTP_SE_nc_wh2_`g': (1 - exp(- (_b[2.wh] / _b[ln_hwage])))), post
		
		* WTP in Euro: children
		est restore coefs_`g'
		eststo WTP_SE_c_`g': nlcom (EUR_WTP_SE_c_fl_`g': (1 - exp(- (_b[fl] + _b[child_fl]) / _b[ln_hwage]))) ///
				(WTP_SE_c_tc_`g': (1 - exp(- ((_b[tc] + _b[child_tc]) / _b[ln_hwage]) ))) ///
				(WTP_SE_c_mw_`g': (1 - exp(- ((_b[mw] + _b[child_mw]) / _b[ln_hwage]) ))) ///
				(WTP_SE_c_wh1_`g': (1 - exp(- ((_b[1.wh] + _b[child_wh1]) / _b[ln_hwage]) ))) ///
				(WTP_SE_c_wh2_`g': (1 - exp(- ((_b[2.wh] + _b[child_wh2]) / _b[ln_hwage]) ))), post
				
				
		** Save WTP in variable (for later: calculation of equivalent wages)
		est restore coefs_`g'
			predictnl wtp_`g'tc = 1 - exp(- (_b[tc]  + _b[child_tc] * child) / _b[ln_hwage]) if woman == `g', ci(wtp_`g'tc_ub wtp_`g'tc_lb)
			predictnl wtp_`g'mw = 1 - exp(- (_b[mw]  + _b[child_mw] * child) / _b[ln_hwage]) if woman == `g', ci(wtp_`g'mw_ub wtp_`g'mw_lb)
			predictnl wtp_`g'fl = 1 - exp(- (_b[fl]  + _b[child_fl] * child) / _b[ln_hwage]) if woman == `g', ci(wtp_`g'fl_ub wtp_`g'fl_lb)
			predictnl wtp_`g'wh1 = 1 - exp(- (_b[1.wh] + _b[child_wh1] * child) / _b[ln_hwage]) if woman == `g', ci(wtp_`g'wh1_ub wtp_`g'wh1_lb)
			predictnl wtp_`g'wh2 = 1 - exp(- (_b[2.wh]  + _b[child_wh2] * child) / _b[ln_hwage]) if woman == `g', ci(wtp_`g'wh2_ub wtp_`g'wh2_lb)

}

** Table: WTP in % and levels by gender
esttab WTP_SE_nc_0 WTP_SE_nc_1 EUR_WTP_SE_nc_0 EUR_WTP_SE_nc_1 using "${dir}\GraphsTables\Table_WTP_Gender_NC.txt", tex varwidth(25) modelwidth(15)  ///
	rename(WTP_SE_nc_fl_1 "Schedule Adaptability" WTP_SE_nc_fl_0 "Schedule Adaptability" EUR_WTP_SE_nc_fl_1 "Schedule Adaptability" EUR_WTP_SE_nc_fl_0 "Schedule Adaptability" WTP_SE_nc_tc_1 Telecommuting WTP_SE_nc_tc_0 Telecommuting EUR_WTP_SE_nc_tc_1 Telecommuting EUR_WTP_SE_nc_tc_0 Telecommuting WTP_SE_nc_mw_1 Meaning WTP_SE_nc_mw_0 Meaning EUR_WTP_SE_nc_mw_1 Meaning EUR_WTP_SE_nc_mw_0 Meaning WTP_SE_nc_wh1_1 "Long Part-Time" WTP_SE_nc_wh1_0 "Long Part-Time" EUR_WTP_SE_nc_wh1_1 "Long Part-Time" EUR_WTP_SE_nc_wh1_0 "Long Part-Time" WTP_SE_nc_wh2_1 "Short Part-Time" WTP_SE_nc_wh2_0 "Short Part-Time" EUR_WTP_SE_nc_wh2_1 "Short Part-Time" EUR_WTP_SE_nc_wh2_0 "Short Part-Time") ///
	nonote ///
    eqlabels(none) collabels(none) mgroups("WtP (\% wage)" "WtP (monthly income, €)", pattern(0 0 1 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span})) mtitles("Men" "Women" "Men" "Women") se nogaps nofloat nonumbers replace b(%5.3f)

esttab WTP_SE_c_0 WTP_SE_c_1 EUR_WTP_SE_c_0 EUR_WTP_SE_c_1 using "${dir}\GraphsTables\Table_WTP_Gender_C.txt", tex varwidth(25) modelwidth(15) b(%5.3f) ///
	rename(WTP_SE_c_fl_1 "Schedule Adaptability" WTP_SE_c_fl_0 "Schedule Adaptability" EUR_WTP_SE_c_fl_1 "Schedule Adaptability" EUR_WTP_SE_c_fl_0 "Schedule Adaptability" WTP_SE_c_tc_1 Telecommuting WTP_SE_c_tc_0 Telecommuting EUR_WTP_SE_c_tc_1 Telecommuting EUR_WTP_SE_c_tc_0 Telecommuting WTP_SE_c_mw_1 Meaning WTP_SE_c_mw_0 Meaning EUR_WTP_SE_c_mw_1 Meaning EUR_WTP_SE_c_mw_0 Meaning WTP_SE_c_wh1_1 "Long Part-Time" WTP_SE_c_wh1_0 "Long Part-Time" EUR_WTP_SE_c_wh1_1 "Long Part-Time" EUR_WTP_SE_c_wh1_0 "Long Part-Time" WTP_SE_c_wh2_1 "Short Part-Time" WTP_SE_c_wh2_0 "Short Part-Time" EUR_WTP_SE_c_wh2_1 "Short Part-Time" EUR_WTP_SE_c_wh2_0 "Short Part-Time") ///
	addnote("Note: standard errors clustered on the individual level and calculated with delta method.") ///
    eqlabels(none) collabels(none) mgroups("WtP (\% wage)" "WtP (monthly income, €)", pattern(0 0 1 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span})) mtitles("Men" "Women" "Men" "Women") se nogaps nofloat nonumbers replace
*/

** T-test: Equality of Coefficients: No Children
forvalues g = 0/1{
	est restore WTP_SE_nc_`g'
	
	scalar wtp_fl_`g' = e(b)[1,1]
	scalar se2_wtp_fl_`g' = e(V)[1,1]
			
	scalar wtp_tc_`g' = e(b)[1,2]
	scalar se2_wtp_tc_`g' = e(V)[2,2]
	
	scalar wtp_mw_`g' = e(b)[1,3]
	scalar se2_wtp_mw_`g' = e(V)[3,3]
		
	scalar wtp_wh1_`g' = e(b)[1,4]
	scalar se2_wtp_wh1_`g' = e(V)[4,4]
	
	scalar wtp_wh2_`g' = e(b)[1,5]
	scalar se2_wtp_wh2_`g' = e(V)[5,5]
}


* P-values
scalar z_fl = (wtp_fl_0 - wtp_fl_1) / sqrt(se2_wtp_fl_0 + se2_wtp_fl_1)
scalar p_fl = normalden(z_fl)

scalar z_tc = (wtp_tc_0 - wtp_tc_1) / sqrt(se2_wtp_tc_0 + se2_wtp_tc_1)
scalar p_tc = normalden(z_tc)

scalar z_mw = (wtp_mw_0 - wtp_mw_1) / sqrt(se2_wtp_mw_0 + se2_wtp_mw_1)
scalar p_mw = normalden(z_mw)

scalar z_wh1 = (wtp_wh1_0 - wtp_wh1_1) / sqrt(se2_wtp_wh1_0 + se2_wtp_wh1_1)
scalar p_wh1 = normalden(z_wh1)

scalar z_wh2 = (wtp_wh2_0 - wtp_wh2_1) / sqrt(se2_wtp_wh2_0 + se2_wtp_wh2_1)
scalar p_wh2 = normalden(z_wh2)

matrix define t_test_nc = (p_fl \ . \ p_tc \ . \ p_mw \ . \ p_wh1 \ . \ p_wh2 \ .)
estadd matrix t_test_nc
esttab matrix(t_test_nc) using "$dir\GraphsTables\P_VAL_t_test_nc.txt", replace

** T-test: Equality of Coefficients: Children
forvalues g = 0/1{
	est restore WTP_SE_c_`g'
	
	scalar wtp_fl_`g' = e(b)[1,1]
	scalar se2_wtp_fl_`g' = e(V)[1,1]
			
	scalar wtp_tc_`g' = e(b)[1,2]
	scalar se2_wtp_tc_`g' = e(V)[2,2]
	
	scalar wtp_mw_`g' = e(b)[1,3]
	scalar se2_wtp_mw_`g' = e(V)[3,3]
		
	scalar wtp_wh1_`g' = e(b)[1,4]
	scalar se2_wtp_wh1_`g' = e(V)[4,4]
	
	scalar wtp_wh2_`g' = e(b)[1,5]
	scalar se2_wtp_wh2_`g' = e(V)[5,5]
}


* P-values
scalar z_fl = (wtp_fl_0 - wtp_fl_1) / sqrt(se2_wtp_fl_0 + se2_wtp_fl_1)
scalar p_fl = normalden(z_fl)

scalar z_tc = (wtp_tc_0 - wtp_tc_1) / sqrt(se2_wtp_tc_0 + se2_wtp_tc_1)
scalar p_tc = normalden(z_tc)

scalar z_mw = (wtp_mw_0 - wtp_mw_1) / sqrt(se2_wtp_mw_0 + se2_wtp_mw_1)
scalar p_mw = normalden(z_mw)

scalar z_wh1 = (wtp_wh1_0 - wtp_wh1_1) / sqrt(se2_wtp_wh1_0 + se2_wtp_wh1_1)
scalar p_wh1 = normalden(z_wh1)

scalar z_wh2 = (wtp_wh2_0 - wtp_wh2_1) / sqrt(se2_wtp_wh2_0 + se2_wtp_wh2_1)
scalar p_wh2 = normalden(z_wh2)

matrix define t_test_c = (p_fl \ . \ p_tc \ . \ p_mw \ . \ p_wh1 \ . \ p_wh2 \ .)
estadd matrix t_test_c
esttab matrix(t_test_c) using "$dir\GraphsTables\P_VAL_t_test_c.txt", replace
*/

** T-test: Equality of Coefficients: Men no Children
foreach c in c nc{
	est restore WTP_SE_`c'_0
	
	scalar wtp_fl_`c' = e(b)[1,1]
	scalar se2_wtp_fl_`c' = e(V)[1,1]
			
	scalar wtp_tc_`c' = e(b)[1,2]
	scalar se2_wtp_tc_`c' = e(V)[2,2]
	
	scalar wtp_mw_`c' = e(b)[1,3]
	scalar se2_wtp_mw_`c' = e(V)[3,3]
		
	scalar wtp_wh1_`c' = e(b)[1,4]
	scalar se2_wtp_wh1_`c' = e(V)[4,4]
	
	scalar wtp_wh2_`c' = e(b)[1,5]
	scalar se2_wtp_wh2_`c' = e(V)[5,5]
}


* P-values
scalar z_tc = (wtp_tc_c - wtp_tc_nc) / sqrt(se2_wtp_tc_c + se2_wtp_tc_nc)
scalar p_tc = normalden(z_tc)

scalar z_mw = (wtp_mw_c - wtp_mw_nc) / sqrt(se2_wtp_mw_c + se2_wtp_mw_nc)
scalar p_mw = normalden(z_mw)

scalar z_fl = (wtp_fl_c - wtp_fl_nc) / sqrt(se2_wtp_fl_c + se2_wtp_fl_nc)
scalar p_fl = normalden(z_fl)

scalar z_wh1 = (wtp_wh1_c - wtp_wh1_nc) / sqrt(se2_wtp_wh1_c + se2_wtp_wh1_nc)
scalar p_wh1 = normalden(z_wh1)

scalar z_wh2 = (wtp_wh2_c - wtp_wh2_nc) / sqrt(se2_wtp_wh2_c + se2_wtp_wh2_nc)
scalar p_wh2 = normalden(z_wh2)

matrix define t_test_m_c_nc0 = (p_tc \ . \ p_mw \ . \ p_fl \ . \ p_wh1 \ . \ p_wh2 \ .)
estadd matrix t_test_m_c_nc0
esttab matrix(t_test_m_c_nc0) using "$dir\GraphsTables\P_VAL_t_test_nc0.txt", replace

** T-test: Equality of Coefficients: Women (no) children
foreach c in c nc{
	est restore WTP_SE_`c'_1
	
	scalar wtp_fl_`c' = e(b)[1,1]
	scalar se2_wtp_fl_`c' = e(V)[1,1]
			
	scalar wtp_tc_`c' = e(b)[1,2]
	scalar se2_wtp_tc_`c' = e(V)[2,2]
	
	scalar wtp_mw_`c' = e(b)[1,3]
	scalar se2_wtp_mw_`c' = e(V)[3,3]
		
	scalar wtp_wh1_`c' = e(b)[1,4]
	scalar se2_wtp_wh1_`c' = e(V)[4,4]
	
	scalar wtp_wh2_`c' = e(b)[1,5]
	scalar se2_wtp_wh2_`c' = e(V)[5,5]
}


* P-values
scalar z_tc = (wtp_tc_c - wtp_tc_nc) / sqrt(se2_wtp_tc_c + se2_wtp_tc_nc)
scalar p_tc = normalden(z_tc)

scalar z_mw = (wtp_mw_c - wtp_mw_nc) / sqrt(se2_wtp_mw_c + se2_wtp_mw_nc)
scalar p_mw = normalden(z_mw)

scalar z_fl = (wtp_fl_c - wtp_fl_nc) / sqrt(se2_wtp_fl_c + se2_wtp_fl_nc)
scalar p_fl = normalden(z_fl)

scalar z_wh1 = (wtp_wh1_c - wtp_wh1_nc) / sqrt(se2_wtp_wh1_c + se2_wtp_wh1_nc)
scalar p_wh1 = normalden(z_wh1)

scalar z_wh2 = (wtp_wh2_c - wtp_wh2_nc) / sqrt(se2_wtp_wh2_c + se2_wtp_wh2_nc)
scalar p_wh2 = normalden(z_wh2)

matrix define t_test_m_c_nc1 = (p_tc \ . \ p_mw \ . \ p_fl \ . \ p_wh1 \ . \ p_wh2 \ .)
estadd matrix t_test_m_c_nc1
esttab matrix(t_test_m_c_nc1) using "$dir\GraphsTables\P_VAL_t_test_nc1.txt", replace


** EXPORT FULL RESULTS (TABLE FOR APPENDIX)
esttab coefs_0 coefs_1 using "${dir}\GraphsTables\OnlineAppendix\Table_Coeffs_Main.txt", tex label ///
	varwidth(25) modelwidth(10) b(%5.3f) noomitted nobase ///
	addnote("Note: standard errors clustered on the individual level.") rename(fl "Schedule Adaptability" tc "Telecommuting" mw "Work Meaning" 1.wh "Long Part-Time" 2.wh "Short Part-Time" ln_hwage "Wages (log)" _cons "Constant")  ///
	eqlabels(none) order("Schedule Adaptability" "Telecommuting" "Work Meaning" "Long Part-Time" "Short Part-Time") collabels(none) mtitles("Men" "Women") se nogaps nofloat nonumbers replace ///
	refcat(child_wh1 "\textbf{Children Interaction}:", nolab)
*/


** EXPORT DATASET WITH PREFERENCES (FOR COUNTERFACTUALS)	
* Grouping gender-specific variables (to later calculate equivalent wages)
foreach i in wh1 wh2 tc fl mw{
	gen wtp_`i' = wtp_0`i' 				if woman == 0
		replace wtp_`i' = wtp_1`i' 		if woman == 1

		}

*** Preparing for export
* One observation per person (back from one per experiment)
bysort nomem_encr: gen obsnr = _n
drop if obsnr != 1

* Generate work hour dummies
gen wh1 = (whcat==1)
gen wh2 = (whcat==2)

* Exporting
save "${dir}\DataCode\Data\Constructed\LISS\Clean_LISS_Preferences.dta", replace
